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ABSTRACT 


A theoretical framework is introduced to model the dynamical changes of the state of polarization during transmission in 
coherent fibre-optic systems. The model generalizes the one-dimensional phase noise random walk to higher dimensions, 
accounting for random polarization drifts, emulating a random walk on the Poincare sphere, which has been successfully 
verified using experimental data. The model is described in the Jones, Stokes and real four-dimensional formalisms, and the 
mapping between them is derived. Such a model will be increasingly important in simulating and optimizing future systems, 
where polarization-multiplexed transmission and sophisticated digital signal processing will be natural parts. The proposed 
polarization drift model is the first of its kind as prior work either models polarization drift as a deterministic process or focuses 
on polarization-mode dispersion in systems where the state of polarization does not affect the receiver performance. We expect 
the model to be useful in a wide-range of photonics applications where stochastic polarization fluctuation is an issue. 


Introduction 

Enabled by digital signal processing, coherent detection allows spectrally efficient communication based on quadrature 
amplitude modulation formats, which carry information in both the intensity and the phase of the optical field, in both 
polarizations. Polarization-multiplexed quadrature phase-shift keying has recently been introduced for 100 Gb/s transmission 
per channel and it is expected that in the near future, higher-order modulation formats will become a necessity to reach even 
higher data rates. However, the demultiplexing of the polarization-multiplexed channels in the receiver requires knowledge 
of the state of polarization (SOP), which is drifting with time. This implies that the SOP must be tracked by a dynamic 
equalizer.*’^ As the SOP changes with time in a random fashion, it is qualitatively different from the chromatic dispersion, 
which can be compensated for using a static equalizer. A deterministic or static behaviour would be straightforward to resolve, 
but with a nondeterministic SOP drift, the equalizer must be continuously adjusted. 

In order to fully understand the impact of an impairment on the performance of a transmission system, a channel model is 
required, which should describe the behaviour of the channel as accurately as possible. Based on the statistical information that 
such models reveal, insights into how to treat the impairments optimally in order to maximize performance can be obtained and 
used as a result. On the other hand, a channel model that does not describe the fibre accurately may hinder the achievement of 
optimal performance. Therefore, it is important that the channel model matches the stochastic nature of the hbre closely. 

Very few results on modelling of random polarization drifts are present in the literature. In the context of equalizer 
development, several models have indeed been suggested, but typically by either using a constant randomly chosen SOP^“® or 
by generating cyclic/quasi-cyclic deterministic changes,^”"* which are usually nonuniform in their coverage of the possible 
SOPs. There is, on the other hand, a wealth of statistical models that describe differential group delay (DGD) and polarization 
mode dispersion dating back to the eighties." However, these results were typically developed to be applicable in systems 
using intensity modulation or single-polarization (differential) phase-shift keying formats, which are not affected by the SOP 
drift. Thus, instead of modelling the time evolution of the SOP, differential equations describing the SOP change with frequency 
and fibre length were typically given. *^’*3 xhere are also some direct measurements of SOP changes, e.g., by Soliman et 
al. However, in their measurement, fast SOP changes were induced in a dispersion-compensating module under laboratory 
conditions, without considering the SOP drift of an entire fibre link. It can be concluded that stochastic SOP drift has so far not 
been given much attention. 

This paper suggests a model for random SOP drifts in the time domain by generalizing the one-dimensional (ID) phase 
noise random walk to a higher dimension. The model is based on a succession of random Jones matrices, where each matrix is 
parameterized by three random variables, chosen from a zero-mean Gaussian distribution with a variance set by a polarization 


linewidth parameter. The latter determines the speed of the drift and depends on the system details. The polarization drift has a 
random walk behaviour, where each step is independent of the previous steps and equally likely in all directions. The model is 
given in the Jones, Stokes formalisms and in the more general real 4D formalism. 

As argued above, the suggested model serves a different purpose than the models of polarization mode dispersion existing in 
the literature. In the latter case, the fibre is viewed either as a concatenation of a small number of segments with relatively high 
DGD, leading to the hinge model,^^’ or the limiting case with segments of infinitesimal lengths, leading to a Maxwellian 
distribution of the DGD.*® The most common assumption is then to model the hinges as independent random SOP rotators.*^’ *^ 
In the anisotropic hinge model, the SOP variation at the hinges is modelled as generalized waveplates parameterized by one 
random angle per hinge.The hybrid hinge model is a further generalized way to model the SOP changes at the hinges,^* but 
in none of these publications is the SOP time evolution discussed. 

The suggested model can be used in simulations for a wide range of photonics applications, where stochastic polarization 
fluctuation is an issue that needs to be modelled. For example, a fibre-optic communication link can be simulated, independently 
of the modulated data as well as of other considered impairments, which can be useful to, e.g., characterize receivers’ 
performance. Moreover, it can reveal statistical knowledge of the received samples affected by polarization rotations, based on 
which the existing tracking algorithms can be optimally tuned or new, more powerful algorithms can be designed. High-precision 
transfer and remote synchronization of microwave and/or radio frequency signals^^ is another application that could benefit 
from a better understanding of how it is affected by polarization drifts, which is currently the limiting factor towards a better 
performance. Other applications such as fibre-optic sensors^^ have been developed for use in a broad range of applications, 
fibre-optic gyroscopes^^ and quantum key distribution^^ are strongly affected by polarization fluctuations and may benefit from 
a better understanding of the transmission medium. 


Present Phase and Polarization Drift Models 

The fibre propagation through a linear medium is often described by a complex 2x2 Jones matrix, which, neglecting the 
nonlinear phenomena, relates the received optical field to the input. Another approach is to use the Stokes-Mueller formalism, 
where the evolution of the Stokes vectors is modelled by a Mueller 3x3 unitary matrix. The latter has the advantage that 
the Stokes vectors are experimentally observable quantities and can be easily visualized as points on a 3D sphere, called 
the Poincare sphere. A further approach to describe the SOP rotations exists in the 4D Euclidean space, *^ where the wave 
propagation can be described by a 4 x 4 real unitary matrix.*® We will focus most of our discussion on the Jones formalism and 
connect it to the Stokes and real 4D formalisms later on. 

An optical signal has two quadratures in two polarizations and can be described by a Jones vector as a function of time t and 
the propagation distance z as E(z,f) = {Ex{z,t),Ey{z,t))'^ , where each element combines the real and imaginary parts of the 
electrical field in the x and y field components and (•)^ denotes transposition. The kth discrete-time input into a transmission 
medium can be written as Ujt = E{Q,IcT), where T is the sampling time and it is commonly related to the symbol interval. The 
received discrete signal, at distance L, is = E{L,kT). 

The propagation of the optical field can be described by a 2 x 2 complex-valued Jones matrix J^, which relates the received 
optical field S C^, in the presence of optical amplifier noise and laser phase noise, to the input S as 

rk = e'^'‘ik\ik + 'ak, ( 1 ) 

where i = s/—l, (j>k models the carrier phase noise and S denotes the additive noise, which is represented by two 
independent complex circular zero-mean Gaussian random variables. Assuming that polarization-dependent losses are 
negligible, the channel matrix Jk can be modelled as a unitary matrix, which preserves the input power during propagation. In 
this work, we assume that the chromatic dispersion has been compensated for and polarization mode dispersion is negligible. The 
transformation J^. can be described using the matrix function J{oik), which is expressed using the matrix exponential [26, p. 165] 
parameterized by three variables a = (ai, a 2 , 0 : 3 ) according to^^’^^ 


J{a) = exp(-/a -d) 

= l 2 cos( 0 ) —/a-CTsin( 0 ), 

where a = (<7i, (T 2 , C 73 ) is a tensor of the Pauli spin matrices 


28 


C7i = 


1 0 
0 -1 


02 = 


0 1 
1 0 /’ 


cr3 = 


0 -i 
i 0 


( 2 ) 


(3) 


This notation of O’,- complies with the definition of the Stokes vector, and it is different from the notation introduced by Frigo.^^ 
The operation a ■ O’ should be interpreted as a scaled linear combination of the three matrices <Ti,< 72 ,<T 3 , and I 2 is the 2x2 
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identity matrix. In general, a 2 x 2 complex unitary matrix has four degrees of freedom (DOFs), but in this case we explicitly 
factored out the phase noise. Including the identity matrix in a would make it possible to account for all four DOFs. The vector 
a can be expressed as a product of two independent random variables a = 0a, i.e., its length 0 = ||a|| in the interval [ 0 , n) and 
the unit vector a = {a\^a 2 ,a^), which represents its direction on the unit sphere. Since Jk is unitary, the inverse can be found by 
the conjugate transpose operation or by negating a^, since = J{—CCk), ■ The same principle holds for the phase noise 

Modern coherent systems require a local oscillator at the receiver in order to get access to both phase and amplitude of the 
electrical held. The local oscillator serves as a reference but it is not synchronized with the transmitter laser, resulting in phase 
noise, which creates the need for carrier synchronization. The phase noise is modelled by the angle (j), while o:i,o :2 and 
model random huctuations of the SOP caused by hbre birefringence and coupling. Both phase noise and SOP drift are 
dynamical processes that change randomly over time. The SOP drift time can vary from microseconds up to seconds, depending 
on the link type. It is usually much longer than the phase drift time, which is in the microsecond range for modern coherent 
systems.^® 

The update rule of the phase noise follows a Wiener process^*'"^^ 

(l>k = ^k + (l>k-i, (4) 

where is the innovation of the phase noise. An alternative form of equation (4) is 

(5) 

which we will generalize later to account for the SOP drift. The innovation (p^ is a random variable drawn independently at 
each time instance k from a zero-mean Gaussian distribution 

(6) 

where Uy = InAvT and Av is the sum of the linewidths of the transmitter and local oscillator lasers. 

The accumulated phase noise at time k is the summation of k Gaussian random variables pi,..., pi; and the initial phase 
po, which becomes a Gaussian-distributed random variable with mean po and variance ka^. The function e‘^'‘ is periodic with 
period 2n, which means that the phase angle pi, can be limited to the interval [0,27r) by applying the modulo 2n operation. In 
this case, the probability density function (pdf) of 0^- becomes a wrapped (around the unit circle) Gaussian distribution. It can 
be straightforwardly verihed that the phase drift model has the following properties; 

1. The innovation pj; is independent of pi for i = 0,... ,k — I. 

2. The innovation is symmetric, in the sense that the probabilities for phase changes p and — (j) are equally likely. 

3. The most likely next phase state is the current state. 

4. The outcome of two consecutive steps e‘^^e“^^ can be emulated by a single step by doubling the variance of p^. 

5. As k increases, the distribution of will approach the uniform distribution on the unit circle. The convergence rate 
towards this distribution increases with the AvT product. 

The time evolution of the pdf of a fixed point corrupted by phase noise is exemplified in Fig. 1. 

The initial phase difference between the two free running lasers has equal probability for every value, therefore it is common 
to model po as a random variable uniformly distributed in the interval [0,271:), i.e., it has equal probability for every possible 
state. 

The autocorrelation function (ACF) quantifies the level of correlation between two samples of a random process taken at 
different time instances by taking the expected value E[ ] of the product of the samples. The ACF of the phase noise with IT 
time separation in response to a constant input u is^^ 


Ar{l) =E[r“ri,+,] 

II 112 / I ^ I 

= l|u|| exp(-—) 

= llu|l^exp(-7r|(|7’Av), 


(7) 
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Figure 1 . Phase noise pdf evolution. The pdf of for 0o = ?r/4 and Gy = 0.0025 is shown. In (a), k—1 corresponds to a 
single innovation and illustrates the second and third properties, i.e., the pdf is symmetric around the current state (the vertical 
line) and the peak of the pdf is at the current state. In (b), k = 5 and the pdf spreads over the circle. In (c), k = 8000 and the pdf 
approaches the uniform pdf, which supports the last property. 


where the operations j-j and H-H denote absolute value and Euclidean norm, respectively. Here we neglected the SOP drift given 
by Ji-, in order to isolate the effects of the phase noise. We will compute the ACE of the SOP drift below. 

In the literature, polarization rotations are generally modelled by the Jones matrix Ji in equation ( 1 ) or subsets of it, obtained 
by considering only one or two of the three DOEs ai,a 2 ,(Xi. Contrary to the phase noise, which is a random walk with respect 
to time, the matrix J^. is in previous literature usually kept constant,^~^’^'^ or it follows a deterministic cyclic/quasi-cyclic rotation 
pattern. The latter is obtained by modelling the parameters of as frequency components (okT. Eor example, aj = COkT 
or 0 ! 2 , eta varied at different frequencies.^’^’ 

Proposed Polarization Drift Model 

We propose to model the polarization drift as a sequence of random matrices, which exploit all three DOEs of J^.. The model 
simulates a random walk on the Poincare sphere and we describe it in the Jones, Stokes and real 4D formalisms. In general, 4D 
unitary matrices have six DOEs, spanning a richer space than the Jones (four DOEs) or Mueller (three DOEs) matrices can. Out 
of the six DOEs, only four are physically realizable for propagating photons, and the remaining two are impossible to describe 
in the Jones or Stokes space.*® The Jones formalism can describe any physically possible phenomenon in the optical fibre 
making it sufficient for wave propagation. The 4D representation is preferred in some digital communication scenarios since 
the performance of a constellation corrupted by additive noise can be directly quantified in this space in contrast to the Stokes 
formalism. Even though the extra two DOEs do not model lightwave propagation, they can be used to account for transmitter 
and/or receiver imperfections, which cannot be done using Jones or Mueller matrices. However, the extra two DOEs are out of 
the scope of this paper and we will focus on the remaining four. 

Jones Space Description 

Similarly to the phase noise update equation (5), we model the time evolution of J^. as 

Ji: =7(aA;)Ji:_i, (8) 

where J{ctk) is a random innovation matrix defined as equation (2). This mathematical formulation is not new to the polarization 
community, as it is commonly used to describe the polarization evolution in space (each Jones matrix represents a waveplate). 
However, here it models the temporal evolution of the SOP, where each innovation matrix corresponds to a time increment. The 
parameters of the innovation J{cCii) are random and drawn independently from a zero-mean Gaussian distribution at each time 
instance k 

otk'^AfiO.a^h), (9) 

where Gp = InApT, and we refer to Ap as the polarization linewidth, which quantifies the speed of the SOP drift, analogous to 
the linewidth describing the phase noise, cf. equation ( 6 ). Drawing a from a zero-mean Gaussian distribution results in special 
cases of 9 and a, where the former becomes a Maxwell-Boltzmann distributed random variable, and the vector a is uniformly 
distributed over the 3D unit sphere, implying that the marginal distribution of each a, is uniform in [—1, 1].^® 
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It is important to note that, contrary to phase noise, the equivalent vector Uk parameterizing in equation ( 8 ) does not 
follow a Wiener process, i.e., <Xk ^ Ctk + Ctk-\, because in general J{ct\)J{ct 2 ) ^ J{ct\ +Ct, 2 )- Equality occurs for (anti-)parallel 
a I and a 2 , and holds approximately when ||Ofi|| ^ 1 and ||a 2 || ^ 1 - In a prestudy for this work,^^ we incorrectly used a 
Wiener process model for the polarization drift. We will return to that model later and discuss its shortcomings. 

Stokes Space Description 

The evolution of the SOP is often analysed in the Stokes space, where the Jones vectors are expressed as real 3D 
vectors and can be easily visualized on the Poincare sphere. The transmitted Jones vector can be expressed as a 
vector [38, eq. (2.5.26)] 

Su^^ufauk (10) 

\E.\^-\Ey\^\ 

2%,{e,e;) , 

-2lm{E,E;)J 

where the ith component of Su^, is given by u^<7,u^ and (•)* denotes conjugation. The equivalent Stokes propagation model of 
equation ( 1 ) can be written 

Sry=MkSuy+Sin. ( 11 ) 

It is important to note that only Sr^, and Suj, can be obtained by applying equation (10) to r^- and VLk, respectively. The noise 
component Sn^ consists of three terms 

Snt = {e'^'^Sk^ikf'Oak + ( 12 ) 

where the hrst two represent the signal-noise interaction and the last one the noise source. It should be noted that there is no 
time averaging in equation (10) and it represents an instantaneous mapping to the Stokes space. Thus, the noise terms are 
polarized and rapidly varying. 

The matrix M<. is a 3 x 3 Mueller matrix, corresponding to the Jones matrix Jk, and the polarization transformation 
introduced by it can be seen as a rotation of the Poincare sphere. It has three parameters, the same as Jk, and can be written as 
Mj. = M{ak), where the function M(-) can be expressed using the matrix exponential^^ 

M{a) — exp(2/C(a)) 

= exp(20/C(a)) (13) 

= I 3 +sin( 20 )/C(a) + (1 - cos( 20 ))/C(a)^ 

where 0 = ||a||, a = a/9 and lC{a) denotes 

( 0 -03 02 \ 

03 0 -01 . (14) 

-02 01 0 / 

The inverse can be obtained by = Mj = M(—a^). The transformation in equation (13) can be viewed in the axis-angle 
rotation description as a rotation around the unit vector a by an angle 20 . 

Note that any operation that applies a phase rotation on both polarizations with the same amount, such as phase noise or 
frequency offsets, cannot be modelled in the Stokes space. This can be seen in equation (11), where the phase noise does not 
alter the transmitted Stokes vector directly, but only contributes to the additive noise in equation (12), which would not exist in 
the absence of n<-. 

Analogously to equation ( 8 ), the time evolution of M^- is 

Mk=M{ak)Mk-u (15) 

where M{ak) is the innovation matrix parameterized by the random vector hk dehned in equation (9). 

Fig. 2 shows an example of an SOP drift as it evolves with time. The line represents the evolution of the vector Mi-(0,0,1)^ 
for k = 1 ,..., 3000. 



Stokes 

Stokes 
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Figure 2. Random walk. The evolution of a random SOP drift obtained by equation (11), without additive noise, for a fixed 
input Su = (0,0,1 and Op =6 - 10^^ is shown. The trajectories for (a) k = 1,..., 300, (b) k = 1, ..., 1500 and (c) 
k— 1,..., 3000 are plotted. 

4D Real Space Description 

Another approach to express the time evolution of the phase and polarization drifts is to use the more uncommon 4D real 
formalism. case, the transmitted/received sample and the noise term in equation ( 1 ) can all be represented as a 

real four-component vector {3ie{Ejc), Im{Ex),3ie{Ey), Im{Ey))^. The transformations induced by the channel are modelled by 
4x4 real unitary matrices. The transformation in equation (1) can be combined into 

Rk ^ R{<pk,0Ck), (16) 

and the function R{-) can be described using the matrix exponential of a linear combination of four basis matrices'® 

R{^,a) = exp(—(/)Ai — ap) 

= exp(—(/)Ai)exp(— a-p) (17) 

= (I 4 cos 0 — A 1 sin 0) (I 4 cos 0 — a • p sin 0), 

where p = {p i,P 2 tP 3 ) and Ai are four constant basis matrices [16, eqs. (20)-(23)] 




-1 

0 

0^ 



0 

0 

-1^ 


f 0 

0 

1 

0^ 


f 0 

1 

0 

o\ 


1 

0 

0 

0 


0 

0 

1 

0 


0 

0 

0 

1 


-1 

0 

0 

0 

Pl = 

0 

0 

0 

1 

; p2 = 

0 

-1 

0 

0 

; P3 = 

-1 

0 

0 

0 

; Ai = 

0 

0 

0 

1 



0 

-1 

V 



0 

0 

0) 


^0 

-1 

0 

V 


^0 

0 

-1 

0; 


(18) 


The inverse can be obtained as = R{—(j)i^, —0,k)- The update of R^. in equation (16) can be expressed analogously to 

equations (8) and (15) as 

^k = R{^k,^k)^k-\, (19) 

where the phase innovation 0^ and the random vector are defined by equation (6) and (9), respectively. 

Polarization Drift Model Properties 

Due to the similarities between the phase noise and the SOP drift, we will use the properties of the phase noise previously 
presented as guidelines to validate the proposed channel model. These properties can be reformulated as follows: 

1. The innovation dfjt is independent of a,- for i = 0,..., k — 1. 

2. The innovation is isotropic, in the sense that all possible orientations of the changes of the Stokes vector corresponding to 
a movement given by one innovation are equally likely. 
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Figure 3. The histograms of M^Su for different steps k and a fixed Su = [1,1,1]^/obtained from the model (top row) and 
from measurements (bottom row) are shown. The highest density is represented by dark red and the lowest by dark blue, the 
outer part of the density. The parameters of the simulated drift in equation (9) are T = 2.2 h (set by the measurements) and 
Ap — 60 nHz (obtained by fitting the dash-dotted ACF line in Fig. 4). In (a) and (d), k — 2 innovation steps are plotted, 
whereas k = 8 in (b), (e) and k = 16 in (c), (f). Gaussian-like isotropic distributions can be noted in all cases, simulations and 
measurements, leading to a good (visual) agreement. The spread over the sphere increases with k and the pdf will become 
uniform if we let k grow large enough. Unfortunately, our measured data do not cover a long enough time period such that 
uniformity is achieved. 


3. The most likely next SOP is the current state. 

4. The outcome of two consecutive steps M{a{)M{a 2 ) can be emulated by a single step M{ai) by doubling the variance 
aj of hk- 

5. As k increases, the pdf of the product M^Su, for a constant Su, will approach the uniform distribution on the Poincare 
sphere. The convergence rate towards this distribution increases with the ApT product. 

In the following, we will investigate these properties and discuss their validity. 

Independent Innovations 

This can be easily concluded from the updating mle in equation ( 8 ), (15) or (19), where the parameters of the innovation do not 
depend on neither the previous innovations nor the state. 

Isotropic Innovation 

We will use the following theorem to show that the innovation M{ct) is isotropic. 

Theorem 1. Let a random unit vector a G /je unifonnly distributed over the 3D sphere, y be a random angle with an 
arbitrary pdf and x gM? an arbitrary unit vector. The pdf of the vector y = M(7a)x, where M{-) is defined in equation ( 13), is 
invariant to rotations around x, i.e., M(j3x)y has the same pdf regardless of 
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In other words, the theorem states that the pdf of the random vector y = M( 7 a)x does not change if it is rotated by any 
angle around x. This can be true only if y is isotropic (centred and equally likely in all directions) around x. The technical 
details of the proof are presented in Supplementary Section I. 

Note that Theorem 1 holds for our proposed Stokes innovation matrix M{cc) since a ~ 9a and a is uniformly distributed 
over the 3D sphere, hence the vector M(d)Su is isotropic around Su- The evolution of M^j-Su can be seen as an isotropic random 
walk on the Poincare sphere starting at MqSu and taking random steps, of size proportional to Up, equally likely in all directions. 
Curiously, the isotropic property is given only by the fact that the rotation axis a is uniformly distributed over the sphere, while 
the angle 9 may have any pdf. In fact, the pdf of the angle 9 determines the shape of the pdf of M(d)Su, which we will discuss 
later. In contrast, our preliminary SOP drift model^^ does not fulhl the isotropicity condition. The update method of the channel 
matrix was done by modelling as Wiener processes, which does not fulhl Theorem 1. 

Pdf of a Random Step 

Unfortunately, we were not able to hnd a closed form expression of the pdf of the point M(d)Su for a hxed Su and a random 
M{a). Using approximations, valid for Up <C 1, which is the case for most practical scenarios, the pdf of M(d)Su can be 
approximated by a bivariate Gaussian pdf centred at Su of variance (T ^\2 on the plane normal to Su. The peak of the pdf is 
located at Su and we can conclude that the next most likely SOP is the current one. The details of the derivations are provided 
in Supplementary Section II. In the same section, under the same assumption of small (7^, we show that the outcome of two 
consecutive random innovations can be achieved by a single random innovation if the variance of the random variables a is 
doubled, which fulhls property 4. 

Uniformity 

The point Mj-Su performs an isotropic random walk on the Poincare sphere, therefore as the number of steps k increases, the 
coverage of the sphere will approach uniformity, meaning that all SOPs will be equally likely. This is intuitive because taking 
random steps in all directions with no preferred orientation will lead to a uniform coverage. This property was observed by 
Zhang et al.,^^ where measurements of a submarine cable resulted in uniform converge of the Poincare sphere. 

We compare the model with experimental results in Fig. 3, where the evolution of the pdf of a Stokes vector affected 
by polarization drift after different numbers of innovations starting from a fixed point is shown. The experimental data was 
obtained by measuring the Jones matrices of a 127 km long buried hbre from 1505 nm to 1565 nm in steps of 0.1 nm for 36 
days at every ^ 2.2 h. The technicalities of the measurement setup and postprocessing have been published elsewhere.^^ The 
histograms corresponding to the measurements were captured from all the Stokes vectors obtained by applying equation (10) to 
the product of the matrix J,(a))JjTi(a)), i.e., the measured matrix corresponding to k innovation steps, with a constant Jones 
vector, where Jr(a)) denotes the measured Jones matrix at time t and wavelength co. 

Initial State 

The initial channel matrix Mq should be chosen such that all the SOPs on the Poincare sphere are uniformly distributed, i.e., 
equally likely, after the Initial step MoSu, for any Su. Analogously, the Initial phase (^o in equation (4) should be chosen from a 
uniform distribution in the interval [0,27r). In order to generate such a matrix Mq, the axis a must be uniformly distributed 
over the unit sphere and the distribution of the angle 9 G [0, ;r/2) must be (1 — cos 9)/7r.^^ The generation of the angle 9 is 
not straightforward, and therefore we present an alternative to generate the axis and the angle simultaneously.'** The vector 
tto = 9a is formed from the unit vector (cos 0,ai sin 0,a2sin0,fl3 sin0)*^ = g/|jg||, where g ^ A/'(0,l4), which will satisfy the 
conditions for both axis a and angle 9. This approach of generating ao of the initial channel matrix can be used regardless of 
the considered method to characterize the polarization evolution, i.e., Jones , Stokes or real 4D formalism. 

The SOP Autocorrelation Functions 

The ACF is commonly used to quantify how quickly, on average, the channel changes In time, frequency, etc.'*^^*' The ACF of 
the SOP drift separated by a time difference of IT in response to a constant input u can be expressed as 

Ar{l) =E[r“ri+/] 

= E[(Jiu)“j^+,u] ^20) 

= |!uf ^(l-(jp)exp(-^)^ , 

where r^. = is the received sample. The details of this derivation can be found in Supplementary Section III. Using the 
approximation (1 — |i|ap/|/|)l*l « exp(—|/|ap), which is valid for (j2< 1 , equation (20) can be approximated as 

A(0«l|ufexp(-3;r|/|rAp). (21) 
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Figure 4. ACF comparison. The normalized ACF of the phase noise and SOP drift is plotted versus normalized time. 
Solid/dashed lines refer to the analytic expressions, whereas the triangles are extracted from measurements. We observe 
excellent agreement between the experiment and theory, except in the tails of the experimental ACF. This inconsistency can be 
caused by the lack of accuracy in that region. 


The ACF expressions in equations (20) and (21) hold for the Jones and real 4D space descriptions. Using an analogous 
derivation as for equation (20), it can be shown that the ACF of the SOP drift in the Stokes space is 

As{l) = 

||2/^2(l-4(72)exp(-2c72) + l 

-IlSull 3 

where Sr^. = M^-Su is the Stokes received sample of a constant input Su affected by polarization drift. Using the following approx¬ 
imations: exp(—2ffp) « 1 — 2(Jp, (Jp « 0 and (1 — S7z\l\TAp/\l\)^'^ « exp(—STrl/lTA/?) for 1, it can be approximated as 

A(/)«||Sufexp(-87r|l|7’Ap). (23) 

Both ACFs of the polarization drift depend only on the time separation I but not on the absolute time k. Comparing equation 
(21) with (23), it is interesting to note that even if they describe the same physical phenomenon, on average, a small movement 
in the Jones/real 4D space will result in a movement that is •\/8/3 larger in the Stokes space. This result was also previously 
observed in similar setups analysing polarization-mode dispersion, where the autocorrelation, with respect to frequency, was 
derived.'*^’^^ The ACF of the phase noise equation (7) has the slowest decay rate, being a factor of three slower than equation 
(21), and a factor of eight slower than equation (23). 

Fig. 4 shows the ACFs of the phase noise and the SOP drift, and the later is compared with experimental results. Here 
the analytic equations (20) and (22) match the experiment for Ap — 80 nHz and Ap — 60 nHz, respectively, and T — 2.2 h. 
We attribute the discrepancy between the two values to the nonunitary of the measured Jones matrices, which, through the 
nonlinear operation in equation (10), cause the inconsistency. 

Linewidth Parameters 

The choice of Av and Ap is important when a system is simulated considering phase noise and/or SOP drifts. Both parameters 
reflect physical properties of the system. The quality of the (transmitter and receiver) lasers can be quantified by the Av 
parameter, which represents the sum of the linewidths of the deployed lasers. Distributed feedback lasers are commonly 
employed in transmission systems due to their cost efficiency. The linewidth of such lasers varies from 100 kHz to 10 MHz.^** 
The polarization linewidth parameter depends on the installation details. Measurements of polarization fluctuations have been 
reported varying from weeks (this work) to seconds,"*^ milliseconds^^’or microseconds under mechanical perturbations.^* 
The polarization linewidth could be quantified through measurements of the ACF, either in the Jones or Stokes space, as in 
Fig. 4. 
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Summary 

This section is intended to summarize the previous sections and provide an easy implementable form of the proposed model 
without requiring knowledge about the details of the derivations. 

First, the parameters of the channel must be selected; 

• The sampling time T, often equal to the symbol interval. 

• The laser linewidth Av, where the contributions of the transmitter and receiver lasers are taken into account. 

• The polarization linewidth Ap. 

The model is then implemented as follows: 

0. Set the initial state of the channel; 

• 0o~W{O,2a:}. 

• Jo = J(®o) using equation (2), where tto = 0a; 0 and a = ( 01 , 02 , 03 ) are identified from the unit vector 
(cos 0,ai sin0,a2sin0,a3 sin0)^ = g/||g||, where g ^ A/’(0,l4). 

1. Fork > 1, update the channel state: 

• + where ~ A/’(0,27rAv7’). 

• Jit = J{cCk)Jk^i, where J(ak) is given by equation (2) and ^ Af (0,271 ApTl 3 ). 

2. Apply the model for every transmitted sample u^, which results in the received sample r^: 

• = e“^'‘JkUk + njc, where denotes the additive noise represented by two independent complex circular zero-mean 
Gaussian random variables. 

Repeat the last two steps for every sample U/t. 

Alternatively, the Jones formalism used above can be replaced by the Stokes formalism using equation (13) instead of 
equation (2) and neglecting the phase noise or 4D formalism by combining and into R*. given by equation (17). In 
the former case, and are 3D Stokes vectors defined as equation (10), whereas in the latter case, f/t n*. are 4D real 
vectors. 

In the summary above, we have included phase and additive noise for completeness. Without loss of generality, the model 
can be applied with/without phase and additive noise; also other impairments can be considered. 

Conclusions 

We have proposed a channel model to emulate random polarization fluctuations based on a sequence of random matrices. The 
model is presented in the three formalisms, Jones, Stokes and real 4D, and it is parameterized by a single variable, called the 
polarization linewidth. 

The model has an isotropic behaviour, which has been successfully verified using experimental data, where every step on 
the Poincare sphere is equally likely in all directions emulating an isotropic random walk and can be easily coupled with any 
other impairments to form a complete channel model. 

Compared to the existing literature, the fundamental advantages of the proposed model are randomness and statistical 
uniformity. Such a model is relevant for a wide range of fibre-optical applications where stochastic polarization fluctuations are 
an issue. It can potentially lead to improved signal processing that accounts optimally for this impairment and more realistic 
simulations can be carried out in order to accurately quantify system performance. 
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I Isotropicity 

We will first prove a lemma, which will be used to prove Theorem 1 in the main paper. 

Lemma 1. For all unit vectors ai, a 2 S and angles 71 , S M we have 

M(7iai)M(72a2) =M(72M(7iai)a2)M(7iai), (SI) 

where M(-) is defined as equation (13). 

Proof. To prove it, we will use the fact that the cross product /C(a)b, where /C(-) is defined in equation (14), is invariant under 


rotations M 

M/C(a)b = ^(Ma)Mb, (S2) 

for all a, b G and any unitary matrix M defined as equation (13). This is true since the vector /C(a)b is orthogonal to a and b 
and its length depends on the area given by the parallelogram formed by a and b. These properties are invariant under rotations, 
hence so is the cross product. Consequently, 

M/C(a) =/C(Ma)M. (S3) 

Now let us denote v = M(7iai)a 2 and, based on equation (S3), we can write 
^(v)M(7iai) =/C(M(7iai)a2)M(7iai) 

= M(7iai)^(a2). (S4) 

By applying equation (S4) twice, we obtain 
^(v)^M(7iai) = ^(v)M(7iai)^(a2) 

= M(7iai)/C(a2)^ (S5) 

The right-hand side of equation (SI) can be simplified as 
M(^M(7iai)a2)M(7iai) =M(72v)M(7iai) 

= (l 3 + sin(272)/C(v) + (l -cos(272))^(v)^)M(7iai) (S6a) 

= M(7iai) + sin(272)^(v)M(7iai) + (l -cos(272))/C(v)^M(7iai) 

= M(7iai) + sin(272)M(7iai )/C(a2) + (1 - cos(272))M(7iai )/C(a2)^ (S6b) 

= M(7iai )(l3 + sin(2^))C(a2) + (1 - cos(272))k:(a2)^) 

= M(7iai)M(72a2), (S6c) 

where equations (S6a) and (S6c) follow from equation (13) and in equation (S6b) we used equations (S4) and (S5). □ 

Now we have the necessary tools to prove Theorem 1 in the main paper. 


Proof. The theorem can be proved by showing that M(j3x)y ^ y for any real angle j3. Let z = M(j3x)y, which can be expressed 
as 

z = M(j3x)M(7a)x 

= M(7M(j3x)a)M(j3x)x (S7a) 

= M(7M(j3x)a)x, (S7b) 

where in equation (S7a) we used Lemma 1. Since the vector a is uniformly distributed over the 3D sphere, the vector M(j3x)a 
is also uniformly distributed over the 3D sphere [1, Def. 6.18], which makes z ~ y. □ 
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II 3D Distribution Anaiysis 

In this section, we derive the approximate pdf of the point Sr = M{cc)Su for a hxed Su and random M{cc). Exact expressions 
are very difficult to obtain, therefore we make use of the approximations sin20 « 20 and cos 20 « 1, valid for (j2 < 1, i.e., 
||a|| ^ 1. Thus M{a) in equation (13) can be approximated as 

M(a) «l 3 + 20/C(a) 

/ 1 —29aj, 20fl2 \ 

« I 20fl3 1 —2001 I 

y—2002 2001 1 / 

/ 1 -2a3 2a2 \ 

« 2a3 1 -2ai . (S8) 

\-2a2 2ai 1 J 


Without loss of generality, we simplify the analysis by setting Su = (1,0,0)^. In this case, based on equation (S8), 
Sr = M(a)Su = ( 1 , 20 : 3 , — 20 : 2 )^ it can be noted that Sr then has a bivariate Gaussian distribution on the plane normal to Su 
and the peak of the distribution centred at Su. 

Using equation (S8) and by removing high order terms, such as a,ay for any i,j, the multiplication of two matrices M{a) 
can be approximated as 


( 1 

M{a)M{^)K 2a3+2)33 

\-2a2-2^2 

KiM{a+P). 


—2(X2 — 2/33 
1 

2o:i + 2j3i 


20:2 + 2)32 \ 

-2ai-2i3i 1 


(S9) 


From equation (S9) we can conclude that two consecutive small innovations can be replaced by a single innovation, i.e., 
M(ai)M(a 2 )Su ~ M(dt)Su, by doubling the variance Cp. 


Ill Autocorrelation 


In this section, we derive the ACF of the SOP drift in equation (20). The derivation uses the Jones description of the model but 
the result is valid for the 4D description as well. 

At hrst we will calculate the expectation of the innovation matrix from equation (2) 

E[7(d)] = E[l2cos(0) — /(fliCTi +a 2 (T 2 + a 3 (T 3 ) sin(0)] 

= E[cos(0)]l2 (SlOa) 

= rcos(0)/e(0)d0l2 
Jo 

= (^(l-(j2)exp(-^)^l2, (SlOb) 

where d = (di,d 2 ,d 3 ) ^ A/’(0,(Jpl 3 ), 0 = ||d|| and a = d/0 = ( 01 , 02 , 03 )- Equation (SlOa) follows because E[o,] = 0 and 
the random variables o;, 0 are independent. Equation (SIOb) follows because the probability density function (pdf) of 0 
is [2, eq. (3.195)] 


fe{0) = 




for 0 > 0. 

The ACF of r*, at time separation I > 0 for a constant input u is 

Arik,k + l)='E[rfrk+i] 

= u^E[J^Ji,+/]u 

= u^E[J^7(d*:+/).. .7(di+i)JA,]u 
= uHE[7(d)]'E[jHj,]u 
= u»E[7(d)]'u. 


(Sll) 


(S12a) 

(S12b) 
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In equation (S12a) we used the fact that the expectation of the innovation matrix is a scaled identity matrix (equation (SlOb)) 
that commutes with J^, and the fact that the innovation matrices J{a,k) are independent. 

Using equation (SlOb) in equation (S12b), the ACF can be expressed as 


Ar{l) 




(SI 3) 


For symmetry reasons, the absolute value of |Z| replaced I in equation (SI 3), making the expression valid for negative / as well. 
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